### alternative RDD specifications ###

# main specification --> kernel weights
rd_mig <- RDestimate(mig ~ date | cntry, data = ESS10, cutpoint = 0, bw = 14)
summary(rd_mig)
rd_eco <- RDestimate(eco ~ date | cntry, data = ESS10, cutpoint = 0, bw = 14)
summary(rd_eco)
rd_cul <- RDestimate(cul ~ date | cntry, data = ESS10, cutpoint = 0, bw = 14)
summary(rd_cul)

lowmig <- rd_mig$ci[1,1]
upmig <- rd_mig$ci[1,2]
coefmig <- rd_mig$est[1]

loweco <- rd_eco$ci[1,1]
upeco <- rd_eco$ci[1,2]
coefeco <- rd_eco$est[1]

lowcul <- rd_cul$ci[1,1]
upcul <- rd_cul$ci[1,2]
coefcul <- rd_cul$est[1]
# slope
slopemig <- lm(mig ~ afg + date + afg*date + cntry, data = subset(smol, wt > 0), weights = wt)
slopeeco <- lm(eco ~ afg + date + afg*date + cntry, data = subset(smol, wt > 0), weights = wt)
slopecul <- lm(cul ~ afg + date + afg*date + cntry, data = subset(smol, wt > 0), weights = wt)

slope1main <- coeftest(slopemig, vcov. = vcovHC(slopemig, type = "HC1"))
slope2main <- coeftest(slopeeco, vcov. = vcovHC(slopeeco, type = "HC1"))
slope3main <- coeftest(slopecul, vcov. = vcovHC(slopecul, type = "HC1"))

slopemigmain <- slope1main[16,1]
upmigmain <- slope1main[16,1] + slope1main[16,2]*1.96
lowmigmain <- slope1main[16,1] - slope1main[16,2]*1.96

slopeecomain <- slope2main[16,1]
upecomain <- slope2main[16,1] + slope2main[16,2]*1.96
lowecomain <- slope2main[16,1] - slope2main[16,2]*1.96

slopeculmain <- slope3main[16,1]
upculmain <- slope3main[16,1] + slope3main[16,2]*1.96
lowculmain <- slope3main[16,1] - slope3main[16,2]*1.96


## manual routine --> population weights

manualw1 <- summary(lm(mig ~ afg + date + afg*date + cntry, data = smol, weights = dweight))
manualw2 <- summary(lm(eco ~ afg + date + afg*date + cntry, data = smol, weights = dweight))
manualw3 <- summary(lm(cul ~ afg + date + afg*date + cntry, data = smol, weights = dweight))

lowmw1 <- manualw1$coefficients[2,1] - manualw1$coefficients[2,2]*1.96
upmw1 <- manualw1$coefficients[2,1] + manualw1$coefficients[2,2]*1.96
coefmw1 <- manualw1$coefficients[2,1]

lowmw2 <- manualw2$coefficients[2,1] - manualw2$coefficients[2,2]*1.96
upmw2 <- manualw2$coefficients[2,1] + manualw2$coefficients[2,2]*1.96
coefmw2 <- manualw2$coefficients[2,1]

lowmw3 <- manualw3$coefficients[2,1] - manualw3$coefficients[2,2]*1.96
upmw3 <- manualw3$coefficients[2,1] + manualw3$coefficients[2,2]*1.96
coefmw3 <- manualw3$coefficients[2,1]
#slope
lowmigmanual <- manualw1$coefficients[16,1] - manualw1$coefficients[16,2]*1.96
upmigmanual <- manualw1$coefficients[16,1] + manualw1$coefficients[16,2]*1.96
slopemigmanual <- manualw1$coefficients[16,1]

lowecomanual <- manualw2$coefficients[16,1] - manualw2$coefficients[16,2]*1.96
upecomanual <- manualw2$coefficients[16,1] + manualw2$coefficients[16,2]*1.96
slopeecomanual <- manualw2$coefficients[16,1]

lowculmanual <- manualw3$coefficients[16,1] - manualw3$coefficients[16,2]*1.96
upculmanual <- manualw3$coefficients[16,1] + manualw3$coefficients[16,2]*1.96
slopeculmanual <- manualw3$coefficients[16,1]

## manual routine --> no weights
manual1 <- summary(lm(mig ~ afg + date + afg*date + cntry, data = smol))
manual2 <- summary(lm(eco ~ afg + date + afg*date + cntry, data = smol))
manual3 <- summary(lm(cul ~ afg + date + afg*date + cntry, data = smol))

lowm1 <- manual1$coefficients[2,1] - manual1$coefficients[2,2]*1.96
upm1 <- manual1$coefficients[2,1] + manual1$coefficients[2,2]*1.96
coefm1 <- manual1$coefficients[2,1]

lowm2 <- manual2$coefficients[2,1] - manual2$coefficients[2,2]*1.96
upm2 <- manual2$coefficients[2,1] + manual2$coefficients[2,2]*1.96
coefm2 <- manual2$coefficients[2,1]

lowm3 <- manual3$coefficients[2,1] - manual3$coefficients[2,2]*1.96
upm3 <- manual3$coefficients[2,1] + manual3$coefficients[2,2]*1.96
coefm3 <- manual3$coefficients[2,1]
#slope
lowmigmanual2 <- manual1$coefficients[16,1] - manual1$coefficients[16,2]*1.96
upmigmanual2 <- manual1$coefficients[16,1] + manual1$coefficients[16,2]*1.96
slopemigmanual2 <- manual1$coefficients[16,1]

lowecomanual2 <- manual2$coefficients[16,1] - manual2$coefficients[16,2]*1.96
upecomanual2 <- manual2$coefficients[16,1] + manual2$coefficients[16,2]*1.96
slopeecomanual2 <- manual2$coefficients[16,1]

lowculmanual2 <- manual3$coefficients[16,1] - manual3$coefficients[16,2]*1.96
upculmanual2 <- manual3$coefficients[16,1] + manual3$coefficients[16,2]*1.96
slopeculmanual2 <- manual3$coefficients[16,1]

## polynomials
tools1 <- rdd_data(ESS10$mig, ESS10$date, cutpoint = 0, covar = ESS10$cntry)
tools2 <- rdd_data(ESS10$eco, ESS10$date, cutpoint = 0, covar = ESS10$cntry)
tools3 <- rdd_data(ESS10$cul, ESS10$date, cutpoint = 0, covar = ESS10$cntry)

# square
square1 <- summary(rdd_reg_lm(tools1, bw = 14, covariates = "ESS10$cntry", order = 2))
square2 <- summary(rdd_reg_lm(tools2, bw = 14, covariates = "ESS10$cntry", order = 2))
square3 <- summary(rdd_reg_lm(tools3, bw = 14, covariates = "ESS10$cntry", order = 2))

lowsq1 <- square1$coefficients[2,1] - square1$coefficients[2,2]*1.96
upsq1 <- square1$coefficients[2,1] + square1$coefficients[2,2]*1.96
coefsq1 <- square1$coefficients[2,1]

lowsq2 <- square2$coefficients[2,1] - square2$coefficients[2,2]*1.96
upsq2 <- square2$coefficients[2,1] + square2$coefficients[2,2]*1.96
coefsq2 <- square2$coefficients[2,1]

lowsq3 <- square3$coefficients[2,1] - square3$coefficients[2,2]*1.96
upsq3 <- square3$coefficients[2,1] + square3$coefficients[2,2]*1.96
coefsq3 <- square3$coefficients[2,1]
# slope
lowsq1slope <- square1$coefficients[5,1] - square1$coefficients[5,2]*1.96
upsq1slope <- square1$coefficients[5,1] + square1$coefficients[5,2]*1.96
coefsq1slope <- square1$coefficients[5,1]

lowsq2slope <- square2$coefficients[5,1] - square2$coefficients[5,2]*1.96
upsq2slope <- square2$coefficients[5,1] + square2$coefficients[5,2]*1.96
coefsq2slope <- square2$coefficients[5,1]

lowsq3slope <- square3$coefficients[5,1] - square3$coefficients[5,2]*1.96
upsq3slope <- square3$coefficients[5,1] + square3$coefficients[5,2]*1.96
coefsq3slope <- square3$coefficients[5,1]


# square double bw
square21 <- summary(rdd_reg_lm(tools1, bw = 28, covariates = "ESS10$cntry", order = 2))
square22 <- summary(rdd_reg_lm(tools2, bw = 28, covariates = "ESS10$cntry", order = 2))
square23 <- summary(rdd_reg_lm(tools3, bw = 28, covariates = "ESS10$cntry", order = 2))

lowsq21 <- square21$coefficients[2,1] - square21$coefficients[2,2]*1.96
upsq21 <- square21$coefficients[2,1] + square21$coefficients[2,2]*1.96
coefsq21 <- square21$coefficients[2,1]

lowsq22 <- square22$coefficients[2,1] - square22$coefficients[2,2]*1.96
upsq22 <- square22$coefficients[2,1] + square22$coefficients[2,2]*1.96
coefsq22 <- square22$coefficients[2,1]

lowsq23 <- square23$coefficients[2,1] - square23$coefficients[2,2]*1.96
upsq23 <- square23$coefficients[2,1] + square23$coefficients[2,2]*1.96
coefsq23 <- square23$coefficients[2,1]
#slope
lowsq21slope <- square21$coefficients[5,1] - square21$coefficients[5,2]*1.96
upsq21slope <- square21$coefficients[5,1] + square21$coefficients[5,2]*1.96
coefsq21slope <- square21$coefficients[5,1]

lowsq22slope <- square22$coefficients[5,1] - square22$coefficients[5,2]*1.96
upsq22slope <- square22$coefficients[5,1] + square22$coefficients[5,2]*1.96
coefsq22slope <- square22$coefficients[5,1]

lowsq23slope <- square23$coefficients[5,1] - square23$coefficients[5,2]*1.96
upsq23slope <- square23$coefficients[5,1] + square23$coefficients[5,2]*1.96
coefsq23slope <- square23$coefficients[5,1]

# cubic
cubic1 <- summary(rdd_reg_lm(tools1, bw = 14, covariates = "ESS10$cntry", order = 3))
cubic2 <- summary(rdd_reg_lm(tools2, bw = 14, covariates = "ESS10$cntry", order = 3))
cubic3 <- summary(rdd_reg_lm(tools3, bw = 14, covariates = "ESS10$cntry", order = 3))

lowcu1 <- cubic1$coefficients[2,1] - cubic1$coefficients[2,2]*1.96
upcu1 <- cubic1$coefficients[2,1] + cubic1$coefficients[2,2]*1.96
coefcu1 <- cubic1$coefficients[2,1]

lowcu2 <- cubic2$coefficients[2,1] - cubic2$coefficients[2,2]*1.96
upcu2 <- cubic2$coefficients[2,1] + cubic2$coefficients[2,2]*1.96
coefcu2 <- cubic2$coefficients[2,1]

lowcu3 <- cubic3$coefficients[2,1] - cubic3$coefficients[2,2]*1.96
upcu3 <- cubic3$coefficients[2,1] + cubic3$coefficients[2,2]*1.96
coefcu3 <- cubic3$coefficients[2,1]
#slope
lowcu1slope <- cubic1$coefficients[6,1] - cubic1$coefficients[6,2]*1.96
upcu1slope <- cubic1$coefficients[6,1] + cubic1$coefficients[6,2]*1.96
coefcu1slope <- cubic1$coefficients[6,1]

lowcu2slope <- cubic2$coefficients[6,1] - cubic2$coefficients[6,2]*1.96
upcu2slope <- cubic2$coefficients[6,1] + cubic2$coefficients[6,2]*1.96
coefcu2slope <- cubic2$coefficients[6,1]

lowcu3slope <- cubic3$coefficients[6,1] - cubic3$coefficients[6,2]*1.96
upcu3slope <- cubic3$coefficients[6,1] + cubic3$coefficients[6,2]*1.96
coefcu3slope <- cubic3$coefficients[6,1]

# cubic double bw
cubic21 <- summary(rdd_reg_lm(tools1, bw = 28, covariates = "ESS10$cntry", order = 3))
cubic22 <- summary(rdd_reg_lm(tools2, bw = 28, covariates = "ESS10$cntry", order = 3))
cubic23 <- summary(rdd_reg_lm(tools3, bw = 28, covariates = "ESS10$cntry", order = 3))

lowcu21 <- cubic21$coefficients[2,1] - cubic21$coefficients[2,2]*1.96
upcu21 <- cubic21$coefficients[2,1] + cubic21$coefficients[2,2]*1.96
coefcu21 <- cubic21$coefficients[2,1]

lowcu22 <- cubic22$coefficients[2,1] - cubic22$coefficients[2,2]*1.96
upcu22 <- cubic22$coefficients[2,1] + cubic22$coefficients[2,2]*1.96
coefcu22 <- cubic22$coefficients[2,1]

lowcu23 <- cubic23$coefficients[2,1] - cubic23$coefficients[2,2]*1.96
upcu23 <- cubic23$coefficients[2,1] + cubic23$coefficients[2,2]*1.96
coefcu23 <- cubic23$coefficients[2,1]
#slope
lowcu21slope <- cubic21$coefficients[6,1] - cubic21$coefficients[6,2]*1.96
upcu21slope <- cubic21$coefficients[6,1] + cubic21$coefficients[6,2]*1.96
coefcu21slope <- cubic21$coefficients[6,1]

lowcu22slope <- cubic22$coefficients[6,1] - cubic22$coefficients[6,2]*1.96
upcu22slope <- cubic22$coefficients[6,1] + cubic22$coefficients[6,2]*1.96
coefcu22slope <- cubic22$coefficients[6,1]

lowcu23slope <- cubic23$coefficients[6,1] - cubic23$coefficients[6,2]*1.96
upcu23slope <- cubic23$coefficients[6,1] + cubic23$coefficients[6,2]*1.96
coefcu23slope <- cubic23$coefficients[6,1]

# main specification --> no FE
nofe1 <- RDestimate(mig ~ date, data = ESS10, cutpoint = 0, bw = 14)
summary(nofe1)
nofe2 <- RDestimate(eco ~ date, data = ESS10, cutpoint = 0, bw = 14)
summary(nofe2)
nofe3 <- RDestimate(cul ~ date, data = ESS10, cutpoint = 0, bw = 14)
summary(nofe3)

lowno1 <- nofe1$ci[1,1]
upno1 <- nofe1$ci[1,2]
coefno1 <- nofe1$est[1]

lowno2 <- nofe2$ci[1,1]
upno2 <- nofe2$ci[1,2]
coefno2 <- nofe2$est[1]

lowno3 <- nofe3$ci[1,1]
upno3 <- nofe3$ci[1,2]
coefno3 <- nofe3$est[1]
#slope
slopemigno <- lm(mig ~ afg + date + afg*date, data = subset(smol, wt > 0), weights = wt)
slopeecono <- lm(eco ~ afg + date + afg*date, data = subset(smol, wt > 0), weights = wt)
slopeculno <- lm(cul ~ afg + date + afg*date, data = subset(smol, wt > 0), weights = wt)

slope1mainno <- coeftest(slopemigno, vcov. = vcovHC(slopemigno, type = "HC1"))
slope2mainno <- coeftest(slopeecono, vcov. = vcovHC(slopeecono, type = "HC1"))
slope3mainno <- coeftest(slopeculno, vcov. = vcovHC(slopeculno, type = "HC1"))

slopemigmainno <- slope1mainno[4,1]
upmigmainno <- slope1mainno[4,1] + slope1mainno[4,2]*1.96
lowmigmainno <- slope1mainno[4,1] - slope1mainno[4,2]*1.96

slopeecomainno <- slope2mainno[4,1]
upecomainno <- slope2mainno[4,1] + slope2mainno[4,2]*1.96
lowecomainno <- slope2mainno[4,1] - slope2mainno[4,2]*1.96

slopeculmainno <- slope3mainno[4,1]
upculmainno <- slope3mainno[4,1] + slope3mainno[4,2]*1.96
lowculmainno <- slope3mainno[4,1] - slope3mainno[4,2]*1.96


# combine coefficients and CIs into dataset

Specification <- c("Kernel Weights (Main)", "Kernel Weights (Main)", "Kernel Weights (Main)", 
                   "Population Weights", "Population Weights", "Population Weights",
                   "No Weights", "No Weights", "No Weights",
                   "2nd Order Polynomial", "2nd Order Polynomial", "2nd Order Polynomial",
                   "2nd Order Double-BW", "2nd Order Double-BW", "2nd Order Double-BW",
                   "3rd Order Polynomial", "3rd Order Polynomial", "3rd Order Polynomial",
                   "3rd Order Double-BW", "3rd Order Double-BW", "3rd Order Double-BW",
                   "Main Without FE", "Main Without FE", "Main Without FE")
Specification <- as.data.frame(Specification)

Outcome <- c("Overall", "Economy", "Culture",
             "Overall", "Economy", "Culture",
             "Overall", "Economy", "Culture",
             "Overall", "Economy", "Culture",
             "Overall", "Economy", "Culture",
             "Overall", "Economy", "Culture",
             "Overall", "Economy", "Culture",
             "Overall", "Economy", "Culture")
Outcome <- as.data.frame(Outcome)

Coefficient <- c(coefmig, coefeco,coefcul,
                 coefmw1, coefmw2, coefmw3,
                 coefm1, coefm2, coefm3,
                 coefsq1, coefsq2, coefsq3,
                 coefsq21, coefsq22, coefsq23,
                 coefcu1, coefcu2, coefcu3,
                 coefcu21, coefcu22, coefcu23,
                 coefno1, coefno2, coefno3)
Coefficient <- as.data.frame(Coefficient)

LowerCI <- c(lowmig, loweco, lowcul,
             lowmw1, lowmw2, lowmw3,
             lowm1, lowm2, lowm3,
             lowsq1, lowsq2, lowsq3,
             lowsq21, lowsq22, lowsq23,
             lowcu1, lowcu2, lowcu3,
             lowcu21, lowcu22, lowcu23,
             lowno1, lowno2, lowno3)
LowerCI <- as.data.frame(LowerCI)

UpperCI <- c(upmig, upeco, upcul,
             upmw1, upmw2, upmw3,
             upm1, upm2, upm3,
             upsq1, upsq2, upsq3,
             upsq21, upsq22, upsq23,
             upcu1, upcu2, upcu3,
             upcu21, upcu22, upcu23,
             upno1, upno2, upno3)
UpperCI <- as.data.frame(UpperCI)

specgraph <- bind_cols(Specification, Outcome, Coefficient, LowerCI, UpperCI)

# plot coefficients and CI

specgraph$Outcome <- factor(specgraph$Outcome)
specgraph$Outcome <- factor(specgraph$Outcome, levels = rev(levels(specgraph$Outcome)))


specgraph1 <- specgraph %>% 
  filter(Specification != "2nd Order Polynomial" & Specification != "2nd Order Double-BW"
          & Specification != "3rd Order Polynomial" & Specification != "3rd Order Double-BW")

levels <- c("Kernel Weights (Main)", "Kernel Weights (Main)", "Kernel Weights (Main)", 
                  "Population Weights", "Population Weights", "Population Weights",
                  "No Weights", "No Weights", "No Weights",
                  "Main Without FE", "Main Without FE", "Main Without FE")

ggplot(specgraph1, aes(x = Specification, y = Coefficient, color = Outcome)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed", size = 1) +
  geom_point(size = 5, position = position_dodge(1)) +
  geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), 
                size = 2, width = 0.4, position = position_dodge(1)) +
  scale_color_manual(values = c("#00BA38", "#00A9FF", "#F8766D")) +
  theme_light(base_size = 35) +
  xlab("Model Specification") +
  ylab("RDD Coefficient") +
  scale_x_discrete(limits = levels) +
  theme(axis.text.x = element_text(angle = 60, hjust= 1))


# plot slopes and CI

specgraph$Slope <- c(slopemigmain, slopeecomain, slopeculmain,
                     slopemigmanual, slopeecomanual, slopeculmanual,
                     slopemigmanual2, slopeecomanual2, slopeculmanual2,
                     coefsq1slope, coefsq2slope, coefsq3slope,
                     coefsq21slope, coefsq22slope, coefsq23slope,
                     coefcu1slope, coefcu2slope, coefcu3slope,
                     coefcu21slope, coefcu22slope, coefcu23slope,
                     slopemigmainno, slopeecomainno, slopeculmainno)

specgraph$Slopeup <- c(upmigmain, upecomain, upculmain,
                       upmigmanual, upecomanual, upculmanual,
                       upmigmanual2, upecomanual2, upculmanual2,
                       upsq1slope, upsq2slope, upsq3slope,
                       upsq21slope, upsq22slope, upsq23slope,
                       upcu1slope, upcu2slope, upcu3slope,
                       upcu21slope, upcu22slope, upcu23slope,
                       upmigmainno, upecomainno, upculmainno)

specgraph$Slopelow <- c(lowmigmain, lowecomain, lowculmain,
                        lowmigmanual, lowecomanual, lowculmanual,
                        lowmigmanual2, lowecomanual2, lowculmanual2,
                        lowsq1slope, lowsq2slope, lowsq3slope,
                        lowsq21slope, lowsq22slope, lowsq23slope,
                        lowcu1slope, lowcu2slope, lowcu3slope,
                        lowcu21slope, lowcu22slope, lowcu23slope,
                        lowmigmainno, lowecomainno, lowculmainno)

specgraph1 <- specgraph %>% 
  filter(Specification != "2nd Order Polynomial" & Specification != "2nd Order Double-BW"
         & Specification != "3rd Order Polynomial" & Specification != "3rd Order Double-BW")

ggplot(specgraph1, aes(x = Specification, y = Slope, color = Outcome)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed", size = 1) +
  geom_point(size = 5, position = position_dodge(1)) +
  geom_errorbar(aes(ymin = Slopelow, ymax = Slopeup), 
                size = 2, width = 0.4, position = position_dodge(1)) +
  scale_color_manual(values = c("#00BA38", "#00A9FF", "#F8766D")) +
  theme_light(base_size = 35) +
  xlab("Model Specification") +
  ylab("Slope Change") +
  scale_x_discrete(limits = levels) +
  theme(axis.text.x = element_text(angle = 60, hjust= 1))






### plot polynomials 

specgraph23 <- specgraph %>% 
  filter(Specification == "2nd Order Polynomial" | Specification == "2nd Order Double-BW"
         | Specification == "3rd Order Polynomial" | Specification == "3rd Order Double-BW")

levels23 <- c("2nd Order Polynomial", "2nd Order Polynomial", "2nd Order Polynomial",
              "2nd Order Double-BW", "2nd Order Double-BW", "2nd Order Double-BW",
              "3rd Order Polynomial", "3rd Order Polynomial", "3rd Order Polynomial",
              "3rd Order Double-BW", "3rd Order Double-BW", "3rd Order Double-BW")

ggplot(specgraph23, aes(x = Specification, y = Coefficient, color = Outcome)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed", size = 1) +
  geom_point(size = 5, position = position_dodge(1)) +
  geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), 
                size = 2, width = 0.4, position = position_dodge(1)) +
  scale_color_manual(values = c("#00BA38", "#00A9FF", "#F8766D")) +
  theme_light(base_size = 35) +
  xlab("Model Specification") +
  ylab("RDD Coefficient") +
  scale_x_discrete(limits = levels23) +
  theme(axis.text.x = element_text(angle = 60, hjust= 1))

ggplot(specgraph23, aes(x = Specification, y = Slope, color = Outcome)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed", size = 1) +
  geom_point(size = 5, position = position_dodge(1)) +
  geom_errorbar(aes(ymin = Slopelow, ymax = Slopeup), 
                size = 2, width = 0.4, position = position_dodge(1)) +
  scale_color_manual(values = c("#00BA38", "#00A9FF", "#F8766D")) +
  theme_light(base_size = 35) +
  xlab("Model Specification") +
  ylab("Slope Change") +
  scale_x_discrete(limits = levels23) +
  theme(axis.text.x = element_text(angle = 60, hjust= 1))



